Phase behavior of a fluid with competing attractive and repulsive interactions 
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Fluids in which the interparticle potential has a hard core, is attractive at moderate separations, 
and repulsive at greater separations are known to exhibit novel phase behavior, including stable 
inhomogeneous phases. Here we report a joint simulation and theoretical study of such a fluid, 
focusing on the relationship between the liquid-vapor transition line and any new phases. The phase 
diagram is studied as a function of the amplitude of the attraction for a certain fixed amplitude 
of the long ranged repulsion. We find that the effect of the repulsion is to substitute the liquid- 
vapor critical point and a portion of the associated liquid-vapor transition line, by two first order 
transitions. One of these transitions separates the vapor from a fluid of spherical liquidlike clusters; 
the other separates the liquid from a fluid of spherical voids. At low temperature, the two transition 
lines intersect one another and a vapor- liquid transition line at a triple point. While most integral 
equation theories are unable to describe the new phase transitions, the Percus Yevick approximation 
does succeed in capturing the vapor-cluster transition, as well as aspects of the structure of the 
cluster fluid, in reasonable agreement with the simulation results. 



I. INTRODUCTION 

Until quite recently, it was commonly supposed that 
the gamut of equilibrium fluid phase behavior exhibited 
by a single component system of particles interacting via 
an isotropic pair potential, does not extend beyond a 
liquid phase and a vapor phase. While this certainly 
appears to be true for prototype models such as the 
Lennard-Jones fluid, it is now recognized that consid- 
erably richer phase behavior can occur [ 1] . For instance, 
systems whose particles have a repulsive core that is suffi- 
ciently 'soft', may exhibit one or more stable liquid-liquid 
phase transitions over and above the liquid-vapor tran- 
sition Furthermore, and notwithstanding their 
isotropy, such soft potentials can give rise to inhomo- 
geneous fluid phases composed of clusters, orpatterned 
morphologies such as stripes and lamellae @, @ S Hi • 

Another class of isotropic pair potential exhibiting rich 
phase behavior is the "mermaid" potential-so called be- 
cause it is attractive at moderate distances and repul- 
sive at large distances (i.e. it has an attractive 'head' 
and a repulsive 'tail') j9|. Effective potentials of the 
mermaid form can be found in colloid-polymer mixtures 
[Tol . [Til . [l2l [l3| : The long ranged repulsion arises from 
the weakly screened charge carried by the colloids, while 
the attraction at intermediate distances stems from the 
depletion forces associated with the non-adsorbing poly- 
mers [l4[. The competition between attraction and re- 
pulsion on distinct length scales is responsible for novel 
behavior such as the appearance of equilibrium cluster 
phases and non-equilibrium gel states [TJ El- Mer- 
maid potentials are further a ppr opriate for the descrip- 
tion of protein solutions [lol Il5j. star-polymer systems 
[lH , the effective interactions between solute particles in 
a subcritical liquid solvent [l7|, and colloidal monolay- 
ers 0, OjjJ]. In the latter (two dimensional) systems, 
the short range attraction arises from the van der Waals 
and capillary forces, while the longer range repulsion is 
thought to be due to dipole-dipole interactions [l8|, [lj| . 




Such two dimensional systems have also been observed 
to form cluster and stripe morphologies [H, [Til H3, HU . 

Beyond their relevance to the real systems listed above, 
models whose particles interact via a mermaid potential 
are of considerable interest from the fundamental per- 
spective of statistical mechanics. Specifically, they elicit 
basic questions such as: What is the scope and charac- 
ter of their phase behaviour, and how can this be de- 
scribed theoretically? In what follows, we briefly review 
the progress to date in addressing these issues. 

A number of theoretical and simulation studies have 
considered aspects of the phase behaviour of a variety of 
fluids interacting via mermaid pot entials [22 

pioneering early work [22|, [23|, [24J , a mean-field (Landau) 
theory approach was developed for systems with compet- 
ing interactions. This predicted that when the amplitude 
of the long ranged repulsion is sufficiently large (relative 
to that of the attraction), modulated phases appear in 
the region of the phase diagram where, if one were only 
to allow for the occurrence of homogeneous phases, the 
theory would predict the liquid- vapor critical point to be 
located. More recent field-theoretic studies [2^, [26|, W\ 
have arrived at similar conclusions. 

Adopting a phenomenological approach, Groenewold 
and Kegel [29|, [30( developed a model to explain how 
competition between short ranged attraction and a longer 
ranged repulsion in colloidal systems could promote clus- 
ter formation. They concluded that depending on the 
relative strengths and ranges of the competing attractive 
and repulsive contributions to the pair potential, large 
clusters (up to several thousands of particles) would be 
stable. Signatures of such clustering were observed by 
Sear and Gelbart within a mean-field (random phase ap- 
proximation (RPA)) liquid state theory 28] study of a 
model in which the attractive and repulsive contributions 
to the potential are both rather long ranged (justifying 
the mean-field approximation). They showed that the 
propensity to clustering is manifest by a peak in the static 
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structure factor, S(k), at a small but non-zero wavevec- 
tor k c . In fact, the mean- field theory predicts a line in the 
phase diagram at which S(k) diverges at k = k c =/= 0, and 
it was inferred that this line indicates microphase separa- 
tion to a modulated phase(s). Such a line has previously 
been dubbed the "A-lin e" [3 2 | in terminology borrowed 
from other contexts [4l|, |42|. |43|| . 

Within the framework of integral equation theory, 
Chen and coworkers [HI, [U, [Hj] solved the Ornstein- 
Zernike (OZ) equation for the structure of a fluid of 
particles interacting via a mermaid potential in which 
the attractive and repulsive contributions are both as- 
signed the Yukawa form (a "double- Yukawa" potential) 
[44|. Two closures were examined: the mean spherical 
approximation (MSA) and the hypernetted-chain (HNC) 
approximation. Focusing on the portion of the phase dia- 
gram where the peak in S(k) at k — k c is developing, the 
authors found that both theories provide a good descrip- 
tion of the fluid structure, as gauged by comparison with 
Monte Carlo (MC) simulations. As will be demonstrated 
in the present work, however, neither the MSA nor the 
HNC is reliable for describing the fluid structure in the 
vicinity of transitions to inhomogeneous phases. 

A generally more accurate integral equation the- 
ory, the self consistent Ornstein-Zernike approximation 
(SCOZA), has been applied in conjunction with hierar- 
chical reference theory (HRT) calculations, by Pini et 
al. [t| [3l[ to a double- Yukawa mermaid potential. The 
authors investigated the influence on the structure and 
phase behaviour of the fluid as the amplitude of the re- 
pulsive contribution was increased from zero. Doing so 
was found to depress the vapor-liquid critical tempera- 
ture and led to the appearance of an anomalously large 
region around the liquid- vapor critical point in which the 
fluid compressibility is very high 0, [3l| . Unfortunately, 
Pini et al. were unable to obtain solutions either from 
the SCOZA or the HRT in the regime where inhomoge- 
neous phases might be expected to occur [H, [3l], H2] • The 
SCOZA results were subsequently compared with a RPA 
density functional theory (DFT) by Archer et al [32j ■ In 
common with earlier mean field approaches [28| , the RPA 
DFT predicts a A-line to occur in the phase diagram when 
the amplitude o f rep ulsive part of the potential exceeds a 
threshold value 32] . In parts of the phase diagram away 
from the A-line, reasonable agreement with SCOZA was 
found for the fluid structure and thermodynamics. 

Sciortino and co-workers [H, H3, HI] have employed a 
variety of theoretical and simulation techniques to inves- 
tigate the properties of mermaid systems. They focused 
attention on potentials having a particularly short ranged 
attractive part (i.e. a small fraction of the particle diam- 
eter) , chosen to mimic the effective pair potential in cer- 
tain colloidal systems. They found - as has been shown 
experimentally [l^] - that the system exhibits both a 
fluid cluster phase and a gel phase comprising intercon- 
nected chains composed of face-sharing tetrahedral clus- 
ters. Indeed, as has been recently emphasized [13, H3], it 
is a general feature of models having particularly short 



ranged attraction, that the microphase separation can oc- 
cur in a regime where the system dynamics are very slow 
and a transition to a non-ergodic state is possible. An- 
other mermaid potential having a similarly short ranged 
attraction was studied using simulation by de Candia et 
al. who observed columnar and lamellar phases [39j . It 
has even been suggested that for such very short ranged 
attractive potentials, the presence of a long ranged re- 
pulsion might not be a prerequisite for cluster formation 

Notwithstanding the extensive body of impressive re- 
sults on a variety of model mermaid systems, displaying 
an intriguing wealth of inhomogeneous phases, central 
questions remain unanswered. Specifically, the detailed 
relationship between the liquid-vapor transition and the 
inhomogeneous phases still seems obscure (recall that 
mean-field theories predict a A-line enclosing a region of 
the phase diagram in which a 'naive' application of the 
theory would predict the liquid vapor critical point to 
be [321].) Furthermore, if the vapor-liquid critical point 
is lost when long ranged repulsion is introduced, what 
happens to the remainder of the liquid-vapor transition 
line? In the present work, we attempt to answer these 
questions by deploying simulation and theory to investi- 
gate the structure and phase behavior of a model system 
whose particles interact via a double- Yukawa mermaid 
potential. 

Our principal findings are as follows. Our Monte Carlo 
simulations show that for a certain (moderately large) 
strength of the repulsive contribution to the pair poten- 
tial, the liquid- vapor critical point is absent. In its stead 
we find two lines of first order phase transitions, each of 
which separates a homogeneous phase from an inhomo- 
geneous (cluster) phase. One of these two transition lines 
is located at low particle number density, and separates 
the vapor from a fluid of spherical liquidlikc clusters; it 
appears to terminate at a critical point at high tempera- 
tures. The other line - located at high density - separates 
a phase of spherical voids from the homogeneous liquid; it 
too appears to terminate at a critical point. At low tem- 
perature, the two transition lines intersect one another 
and a vapor-liquid transition line at a triple point. 

We complement our simulation studies with an investi- 
gation of the utility of a number of standard liquid state 
theories for describing the phase behaviour of our model. 
We first apply the DFT of Ref. [HJ to trace the locus of 
the A-line, noting that if one interprets this line as repre- 
senting the transition to periodically modulated phases, 
then its topology is incompatible with the phase diagram 
emerging from the simulations. Turning our attention to 
integral equation theories, we find that the HNC has a no- 
solution region in the portion of the phase diagram where 
the transitions to inhomogeneous phases occur, although 
it does yield a portion of the vapor-liquid transition. Use 
of a simple modified HNC (MHNC) approximation sim- 
ilarly fails to provide a solution in the region of interest. 
We further find that (for separate reasons) the MSA is 
of little use in this region of the phase diagram. Inter- 
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estingly, however, the Percus Yevick (PY) approximation 
is able to describe the vapor-cluster phase transition as 
well as key aspects of the structure of the cluster fluid, 
in reasonable agreement with the MC results. 

Our paper is organized as follows. We introduce our 
model mermaid potential in Sec. |TTJ The MC simulation 
methodology used to study the model, and our findings 
concerning the phase behaviour and the character of the 
inhomogeneous phases are described in Sec. IHII An in- 
vestigation of the utility of mean field DFT and integral 
equations for describing the phase behavior is detailed 
in Sec. IIVI Finally, a discussion of the implications of 
our findings and the outlook for future work features in 
Sec. El 



II. MODEL 

The model that we have elected to study comprises 
an isotropic two-body interparticle potential of the hard- 
core plus double- Yukawa form: 

{oo r < a 

-ecrexp(-zi(r/cr — l))/r (1) 
+A<jexp(— z 2 (r/a — l))/r r > a. 

Here f3 = l/fc^T is the inverse temperature, which we set 
equal to unity, while a is the particle diameter. The first 
Yukawa term represents an interparticle attraction whose 
strength is controlled by a parameter e > 0, while the 
second term, with strength parameter A > 0, represents 
a repulsion. We shall focus on the regime in which the 
range of the repulsion exceeds that of the attraction i.e. 
z\ > z 2 , leading to a repulsive tail to the potential. 

The specific choice of z\ and z 2 requires a balance 
to be struck between a number of competing desider- 
ata. Firstly, the attractive part should not be so short- 
ranged that the equilibrium liquid-vapor phase behavior 
is wholly preempted by freezing |46[ and/or dynamic 
arrest (gelation) (4p| , phenomena which - while undoubt- 
edly of interest in their own right - are not the specific 
concern of this work. Secondly, the range of the repulsive 
part must exceed that of the attractive contribution, and 
be sufficiently large to enable meaningful comparisons 
with mean field theories, the accuracy of which increases 
with the potential range. Finally, from a purely practical 
standpoint, the range should not be so great that very 
large system sizes are necessary in order to obviate cutoff 
artifacts when the potential is truncated in the interests 
of computational efficiency. A satisfactory compromise in 
these respects was found by assigning Z\ = 2, z% = 1, and 
truncating the potential at r = r c = 7.0a. For reasons 
relating to the occurrence of large length-scale inhomo- 
geneous phases, no mean field correction was applied for 
the effects of the potential truncation. 

As described in the Introduction, physically one can re- 
gard the potential of Eq. (fTJ as capturing the combined 



overall effect of a screened coulombic repulsion between 
charged colloidal particles, and a shorter-ranged attrac- 
tion engendered by a depletion agent such as polymer 
chains. In this spirit, it is natural to regard e as a measure 
of the polymer fugacity, A as a measure of the colloidal 
charge, and treat both parameters independently. In the 
present work, however, we shall primarily be concerned 
with the phase behavior as a function of e _1 at fixed A. 
For the latter, we have assigned the value A = 0.55. This 
choice was motivated by a preliminary study of the model 
within DFT [33] , which suggested that (given the choice 
Z\ = 2,z 2 = 1) this strength of repulsion is sufficient to 
engender a A-line encompassing the critical region of the 
liquid-vapor line. One might therefore hope that qualita- 
tive alterations to the standard scenario of liquid-vapor 
phase behavior would ensue. 

The form of the potential f3v{r) is shown in Fig. [T] for 
parameter values e^ 1 = 0.4, A = 0.55, Z\ = 2, z 2 = I. 
Also shown is the form of (r /a) 2 f3v(r), which provides 
a useful indication of which cutoff values, r c , are likely 
to result in significant corrections to the internal energy 
compared to the full potential. The figure confirms that 
for our choice of the cutoff, r c = 7. Oct, only small cor- 
rections to the limiting behaviour are to be expected. It 
should be noted, however, that because this cutoff value 
is much larger than the values typically used in the sim- 
ulation of simple fluids such as Lennard-Jonesium, any 
simulation study is expected to entail a considerable com- 
putational investment. 
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FIG. 1: The form of the potential f3v(r) (Eq. [TJ for param- 
eters e _1 = 0.4,^4 = 0.55,2i = 2,22 = 1. Also shown for 
comparison is the form of (r /a) 2 (3v(r). 



III. SIMULATION STUDIES 
A. Techniques 

We have studied the phase behaviour of the model 
of Sec. |TT] using a grand canonical MC simulation algo- 
rithm [47| . Where possible, accurate location of points 
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of first order phase coexistence was facilitated by the use 
of the multicanonical preweighting technique [48| , aided 
by multihistogram reweighting [49| according to the pro- 
cedure described in ref . [50( . Most of the results we shall 
present were obtained for a system of linear size L — 21er, 
although some data was also collected for L = 28a in 
order to gauge the scale of finite-size effects. Periodic 
boundaries were employed throughout. 

As shall be described below, techniques were imple- 
mented to determine the distribution of sizes of particle 
clusters in selected regions of the phase diagram. A clus- 
ter comprises a subset of particles that are interlinked via 
pathways of interparticle 'bonds'. However, in contrast 
to lattice models, the definition of a bond in a system 
with continuous translational symmetry is somewhat am- 
biguous. We adopt a criterion which derives from that 
used for cluster identification in spin models. Specif- 
ically, we determine the interaction energy v between 
each pair of particles and assign a bond with probability 
Pbond = 1 — exp(/3v). Clusters of bonded particles are 
then identified using the efficient enumeration algorithm 
of Hoshen and Kopelman . 



B. Phase diagram 

Prior to commencing exploration of the phase diagram, 
an initial estimate was required for the range of values 
of the attractive strength for which non-homogeneous 
states might be expected to occur. Guided by the DFT 
calculations reported in Sec. IIVI (which predict that for 
A = 0.55, a A-line appears for 0.54 > e _1 > 0.48), we 
selected e^ 1 — 0.5 as a suitable candidate. Subsequently 
it was found that smaller values of e _1 were necessary to 
generate inhomogeneous phases. 

For each value of e _1 studied, the dependence of the av- 
erage of the fluctuating particle number density, p, on the 
applied chemical potential p was measured. The result- 
ing forms of p(p) are shown in Fig.[2fa). From this figure 
one observes that for the initial value e _1 = 0.5, p{p) is 
essentially smooth, having a form reminiscent of a one- 
phase (super-critical) fluid. However, as e _1 is reduced 
from this value, the gradient of the curves for moder- 
ate densities gradually increase in magnitude, indicating 
an increase in the fluid compressibility. Concomitantly, 
kinks in the curves start to develop at low and high den- 
sity, which sharpen into discontinuities as e _1 is further 
decreased. For e _1 < 0.42, these low and high density 
kinks are supplemented by additional ones at intermedi- 
ate densities, as is apparent from the close-up of the data 
for e -1 = 0.41 and e~ 1 = 0.4 shown in Fig.^b). 

We address first the phenomena underlying the ap- 
pearance of the discontinuities in p(p) at low and high 
density. These arise from first order phase transitions. 
In the low density case, the transition is between a vapor 
and an inhomogeneous phase composed of large spherical 
liquidlike clusters; in the high density case it is from a 
liquid to an inhomogeneous liquid phase containing large 
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FIG. 2: (a) Measured forms of p(p) for zi = 2, Z2 = 1, 
A — 0.55 and various e _1 . (b) An enlargement showing the 
data for e _1 = 0.41 and e~ 1 = 0.4 as described in the text. 



spherical vaporlike voids (i.e. bubbles). Discussion of the 
character of these inhomogeneous structures is deferred 
until Sec. IhTCI 

Evidence for the existence of these phase transitions 
comes from the measured forms of the distribution of 
the fluctuating instantaneous number density p(p) at the 
model parameters for which the kinks occur. For the 
low density transition, these distributions are shown in 
Fig. [3ja) . One observes a two-peaked structure; the nar- 
row peak at low densities corresponds to the vapor phase, 
while a much broader higher density peak corresponds to 
the cluster phase. At small values of e" 1 the two peaks 
are widely separated, the trough between them is deep 
(cf. the log scale of fig-[3Jb)), and the cluster peak is very 
broad (this latter feature reflects the large compressibil- 
ity associated with the steepness of p(p) at these values 
of -recall Fig. [D(a)). As e _1 is increased, the two 
peaks approach one another and eventually merge into a 
single peak; this occurs at a value of e _1 consistent with 
that at which the sharp low density kink in p{p) becomes 
smoothed out. 

The computational cost of performing long simulations 
in the region of the high density transition, verges on the 
prohibitive. Consequently, we have been able to obtain 
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FIG. 3: (a) The distribution of the fluctuating number den- 
sity for state points corresponding to the low density vapor- 
spherical cluster transition, (b) The same data plotted on a 
log scale. 
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FIG. 4: The distribution of the fluctuating number density 
for the state point e" 1 = 0.41, f3p, = —3.91, lying on the high 
density spherical bubble- liquid transition. 
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the coexistence form of p(p) for only a single value of e , 
namely e -1 = 0.41. This distribution is shown in Fig. [3] 
and exhibits a narrow liquidlike peak at high density to- 
gether with a much broader peak at lower densities which 
corresponds to a liquid containing a large void, as will be 
discussed in Sec. IIII CI Again the parameters at which 
we find the double-peaked structure match those of the 
corresponding kink in p(p) , which can therefore be taken 
as an alternative signature of the transition. 

Unfortunately, it proved impossible to accurately 
quantify the number and energy densities of the clus- 
ter phases that coexist with the respective vapor and 
liquid phases. The problems are traceable to the large 
length scale associated with the typical cluster size, as 
illustrated in Fig. [5fa) , which shows the measured den- 
sity distribution p{p) for two system sizes at a near- 
coexistence point located well inside the two-phase vapor- 
cluster region. From the figure, it is evident that while 
the position of the vapor phase peak is essentially system 
size- independent, the density of the cluster peak shifts 
strongly to lower density as L is increased. Thus the sys- 
tem sizes attainable in this work fail to provide access 



FIG. 5: Near-coexistence density distributions for (a) e" 1 — 
0.415 and (b) e _1 = 0.425. In each case data is shown for two 
system sizes. 



to the thermodynamic limit for this phase. Analogous 
effects are found for the energy distribution (data not 
shown) . 

The low and high density first order transition lines in 
e — fx space are presented in Fig. [BJ In this figure, the 
estimates of points on the low density transition line de- 
rive from the data of Fig. [3] by locating the values of fi 
for which p{p) exhibits two peaks of approximately equal 
weight. Those for the high density transition derive both 
from the data of Fig. |3] for the case e _1 = 0.41, and oth- 
erwise from the position of the high density kink in the 
p(p) curves. The loci of the two transition lines delimits 
a region within which inhomogeneous structures occur. 
For small e" 1 ^ 0.39, this region tapers down to a point, 
beyond which we were unable to stabilize inhomogeneous 
phases, instead finding a single transition from a vapor 
to a liquid. This suggests that a triple point occurs at 
e" 1 ^ 0.39 below which standard vapor-liquid coexistence 
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FIG. 6: Estimates of the phase diagram, as described in 
the text. Circles are the vapor-spherical cluster transition, 
squares are the spherical bubble- liquid transition. Crosses lie 
on the liquid-vapor coexistence line, lines merely guide the 
eye. Filled symbols locate the putative critical points. Un- 
certainties do not exceed the symbol sizes. 



occurs. We have estimated the locus of a portion of the 
vapor-liquid coexistence from the center of the hysteresis 
loop in e — ji formed by traversing the transition from 
vapor to liquid and back again, and these estimates are 
also marked on Fig. [5] Owing to the great number of 
particles in the liquid phase, the considerable potential 
cutoff distance, and the substantial density difference be- 
tween the liquid and vapor, it was not possible to link the 
phase spaces of coexisting vapor and liquid directly us- 
ing biased sampling techniques, as was done for the low 
density (vapor-cluster) transition. 

As regards the nature of the transitions for large e , 
one observes from Fig. [3] that for the vapor-cluster tran- 
sition, the two peaks in p(p) coalesce at a value of 
close to that at which the low density kink in p(p) dis- 
appears. It seems reasonable to assume that a critical 
point occurs in this region and that, by extension, the 
high density transition similarly end at a critical point 
near to where the high density kink in p(p) vanishes. 
However, accurately pinning down the respective critical 
point parameters, proved problematic. The source of the 
difficulty is, again, the large cluster length scale. Specifi- 
cally, on increasing e _1 , the two peaks merge before any 
system size dependence of the vapor peak becomes ap- 
parent. This implies that (in contrast to the case for e.g. 
a liquid- vapor transition in a simple fluid (H^ j) the double 
peaked structure of p(p) is lost before the thermal corre- 
lation length becomes comparable with the system size. 
Accordingly, the value of e _1 at which the peaks merge 
in a finite-sized system constitutes an underestimate of 
the true critical point value of e _1 and can thus, strictly 
speaking, only provide a lower bound on the critical point 
value of e _1 [53|. The existence of an additional large 
length scale for one of the phases severely complicates 
the implementation of standard finite-size methods for 



locating criticality, which are based on the assumption of 
a single dominant length scale in the near-critical region, 
namely the thermal correlation length. Nevertheless, for 
e _1 > 0.43 the low and high density kinks in p(p) have 
disappeared, the curves are smooth, and we thus tenta- 
tively (and conservatively) assign 0.425 < e^ 1 < 0.440 
for both the low density and high density transitions. 



C. Nature of the coexisting phases at the low and 
high density transitions 

In Sec. lIII 51 we mapped the boundary between the ho- 
mogeneous and inhomogeneous phases. Here we ponder 
the character of the inhomogeneous phases in more detail. 
Fig. EJa) and (b) show, for e _1 = 0.41, typical configura- 
tion snapshots of the inhomogeneous phase which coex- 
ists with homogeneous vapor (liquid) at the low (high) 
density transitions respectively. For the low density tran- 
sition, a vapor coexists with a phase containing an ap- 
proximately spherical cluster, whose local density is liq- 
uidlike. For the high density transition, a dense liquid 
coexists with a phase containing a large spherical bubble 
of vapor. 

The visual identification of clusters at the low density 
transition is corroborated by measurements of the radial 
distribution function g(r) of the coexisting phases. Fig. [5] 
displays the measured forms of g(r) for the vapor and 
the coexisting spherical cluster phase at e _1 = 0.41. One 
notes that while the vapor phase is, to a great extent, 
structureless (i.e. g(r) effectively reaches its asymptotic 
value of unity for r > 3c), the cluster phase exhibits un- 
usual features, namely a very pronounced enhancement 
in the value of g(r) extending over a considerable length 
scale that is indicative of the cluster radius. Furthermore, 
as the inset shows, g(r = 10.5er) ~ 0.8 -significantly 
less than unity. This latter feature heralds the onset 
of large- length-scale oscillations in g(r) which have their 
origin in the inter-cluster correlations of the cluster fluid 
phase [32|]. Further evidence for this assertion derives 
from our integral equation calculations to be presented 
in Sec. HVBl 

The question naturally arises as to how the character of 
the spherical cluster phase alters as one tracks the transi- 
tion line by varying e _1 . To answer it, in part at least, we 
have implemented a cluster identification algorithm [Hit ] 
and used it to obtain the form of the cluster mass distri- 
bution P(s) at the vapor-cluster transition, for a number 
of coexistence state points. The results (Fig. ^ show, in 
each instance, a double peaked distribution. The peak 
at low mass occurs at s = 1, indicating that little or no 
clustering occurs in the vapor. The broad peak occurring 
at high mass indicates that the liquidlike clusters contain 
several hundred particle. The latter peak shifts strongly 
to greater masses with decreasing e~ 1 , showing that the 
typical radius of the clusters grows accordingly; this fea- 
ture has important implications for finite-size effects, as 
we discuss below. 
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FIG. 9: Cluster mass distribution p(s) (shown on a log scale) 
for a selection of points along the vapor-cluster transition. 
The curves have been shifted vertically to aid their distin- 
guishability. 




D. Further inhomogeneous states 



FIG. 7: (a) A typical configurations of the spherical cluster 
phase which coexists with vapor (not shown) at the low den- 
sity transition, (b) A typical configuration of the spherical 
vapor bubble phase that coexists with homogeneous liquid 
(not shown) at the high density transition. 



Finally in this section, we record that we have not 
attempted to characterize in any detail the properties 
of the inhomogeneous vapor-bubble phase occurring at 
the high density transition, due to the prohibitively large 
computational effort that this would entail. 



We consider next the structures that form as one tra- 
verses the region of inhomogeneous states separating the 
low and high density transitions. On increasing /i at con- 
stant e~ 1 = 0.41 from its coexistence value at the low den- 
sity (vapor-cluster) transition, visual inspection of typical 
configurations (Fig. [T0|) shows that the spherical clusters 
are replaced firstly by cylindrical clusters (which span the 
periodic system in one direction), and then, at still higher 
densities (per 3 ~ 0.2), by slablike structures (which span 
in two directions). In our simulations, one structure ap- 
peared to evolve smoothly into the next, and we could 
discern no clear signature of discontinuities in the gradi- 
ent of the p(fi) in this range of densities (cf. Fig. (21(b)). 
Hence on this basis, there is no evidence for the exis- 
tence of first order phase transitions between spheres and 
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cylinders or cylinders and slabs. However, given that the 
characteristic length scales of these structures is compa- 
rable with our system size, we cannot rule out that sharp 
transitions could become apparent for much larger sys- 
tem volumes. 

Moving to higher densities (pa 3 > 0.3), two sharp 
kinks are visible in the p(p) curve for (T 1 — 0.41 (see 
Fig. [DJb)). The first kink, appearing at a density of 
pa 3 ~ 0.37, occurs in a region in which the configurations 
are slab-like, and appears to relate to a single cuboidal 
slab being replaced by two parallel slabs. Accommodat- 
ing an additional slab within the simulation box boosts 
the number of particles in each that are within range of 
the long ranged repulsions from particles occupying the 
other slab. Consequently it becomes energetically less 
favorable for the slabs to thicken as p is increased fur- 
ther, leading to a decrease in the compressibility (and 
thus too, the gradient of p(p))- Accordingly, the kink at 
pa 3 ~ 0.37 is something of a finite-size artifact. Visual- 
ization of configurations either side of the second inter- 
mediate density kink (occurring at pa 3 ~ 0.47) reveals 
that this is associated with a transition from the parallel 
slabs of liquid just described, to a liquid containing large 
cylindrically shaped voids, or "bubbles" of vapor which 
span the system in one direction (Fig. HOf c)). A jump in 
the density at pa 3 ~ 0.7 corresponds to cylindrical voids 
being replaced by spherical voids. 

Turning now to the data for e _1 = 0.4, shown in 
Fig. [2Kb) , the kinks that form the plateau in p(p) at mod- 
erate densities have an origin distinct from that just de- 
scribed for the case e _1 = 0.41. The configurations either 
side of the first kink contain a single liquidlike slab. The 
plateau seems to develop because -as the density grows- 
the two surfaces of this single slab are forced together via 
the periodic boundaries until, at some point, the particles 
near one surface come within the range of the repulsive 
part of the potential of the particles at the other surface. 
Thus again, this feature would appear to be a finite-size 
artifact. The second kink marks an abrupt condensation 
of the slab into the homogeneous liquid phase. 

The character of the inhomogeneous structures that 
occupy the region separating the pure phases helps to 
explain the dramatic increase in the compressibility (and 
hence the gradient of p(p)) on reduction of e . This 
increase is (we believe) largely attributable to the in- 
creasing prevalence of slab like structures as e _1 is re- 
duced. Since the size of the spherical clusters that coex- 
ist with the pure phases grows with decreasing e — x , they 
eventually approach that of the system size, at which 
point spheres or cylinders are no longer stable and are 
replaced by slabs. Since it costs little free energy to move 
a slab interface (except when the two interfaces interact 
via the periodic boundaries), the compressibility grows 
very large. We further note that for values of e _1 < 0.39, 
i.e. below that of the triple point, one should expect such 
behavior since at densities intermediate between those of 
the stable coexisting vapor and liquid phases, slab-like 
configurations dominate [54j . 




FIG. 10: Snapshots of the typical sequence of finite system 
size-limited structures that were observed for e _1 = 0.41 when 
traversing the region of densities intermediate between the 
density of the spherical liquidlike cluster phase (which coexists 
with the vapor cf. Fig-Efa)), and the density of the spherical 
void phase (which coexists with the liquid cf. Fig. [7{b)). (a) 
A cylinder (b) A slab (c) A cylindrical void. See text for 
more details. 

Summing up the findings of this subsection, we have 
presented evidence for a variety of inhomogeneous struc- 
tures: spherical and cylindrical liquidlikc clusters, sin- 
gle and multiple liquidlike slabs, cylindrical and spherical 
bubbles. However, the precise sequence of structures that 
occurs when traversing the inhomogeneous region at fixed 
e _1 depends both on the value of (T 1 itself, and on the 
system size. This is perhaps not altogether surprising, 
given that the typical length scale of the inhomogeneous 



structures depends on (T 1 (cf. Fig. [9]), and that for suffi- 
ciently small e _1 , this length scale can exceed the linear 
extent of the simulation box. Additionally, it should be 
mentioned that protracted relaxation times were encoun- 
tered in the inhomogeneous region of the phase diagram. 
These arise because a local (i.e. single particle) update 
procedure requires many iterations in order to decorre- 
late a system having an inherently large length scale. Ac- 
cordingly, great computational expenditure is necessary 
to ensure that the system attains the equilibrium struc- 
ture for a given state point. In fact, we observed a certain 
amount of irreproducibility regarding the precise form of 
p(pr) obtained on traversing the inhomogeneous region 
from vapor to liquid compared to the reverse path. The 
source of this irreproducibility is presumably traceable, 
in part at least, to the extended relaxation times. 
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IV. MEAN FIELD AND INTEGRAL EQUATION 
THEORETICAL STUDIES 

A. Mean-field theory of the fluid structure and 
thermodynamics 

In Ref. [32| the authors developed a mean-field DFT 
theory within the random phase approximation (RPA), 
for systems interacting via potentials of the form in Eq. 

We will not describe in detail the theory here - 
instead referring interested readers to Ref. [32j ■ The key 
idea behind this approach (and indeed most other mean- 
field approaches) is to split the pair potential into two 
contributions: a reference part, v r (r) and the remainder 
or 'perturbation', v p (r), i.e. v(r) — v r (r) + v p (r) . For the 
present system an obvious choice for the reference part 
is the hard sphere potential, v r (r) = Wfc s (r), where 
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and therefore: 



r < a 
r > a, 



(2) 



FIG. 11: Phase diagram for zi = 2, zi = 1 and A = 0.55, 
obtained from the RPA DFT theory of Ref. [H . The solid line 
is the binodal and the dashed line is the A-line. (a) Plotted 
in the density p versus e _1 plane, (b) plotted in the chemical 
potential /i versus e _1 plane. 
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r < g 
r > a. 



(3) 



Note that there is no unique choice for v p (r) in the range 
r < a, inside the hard core. The choice made in Eq. (J3j> 
is the same as that of Ref. [32[ . 

The Helmholtz free energy of the system is approxi- 
mated as follows 1321: 



F\p]~F«r\p] + ± J&vJ drV(r)p(r> p (|r-r'|), (4) 

where p(r) is the one body density profile of the fluid and 
Fhs° s [p] i s t ne Rosenfeld approximation for the Helmholtz 
free energy of a hard sphere fluid [H, [56|, [57| • For the bulk 



fluid this amounts to the following (RPA) approximation 
for the pair direct correlation function [321 . |44| : 



c(r;p) 



■ Pv p (r), 



(5) 



where c^ s (r; p) is the Percus-Yevick approximation for 
the hard-sphere direct pair correlation function for a bulk 
fluid of density p [321 . |44| . 

For the potential parameters of concern in the present 
work (z\ = 2, Z2 = 1 and A = 0.55), the mean-field the- 
ory of Ref. [33] yields the phase diagram displayed in Fig. 
[TTI From this figure one sees that the theory predicts a 
A-line enclosing the region of the phase diagram contain- 
ing the liquid- vapor critical point. The A-line is defined 
as the locus of points in the phase diagram for which 
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the static structure factor S(k) diverges at a particular 
wave number k c ^ [3^|. We take the A-line to indicate 
that the theory predicts a phase transition to a periodi- 
cally modulated inhomogeneous phase. This phase pre- 
empts the standard liquid vapor behavior that normally 
occurs in the subcritical region of a simple fluid. We 
note, however, that the topology of the A-line and hence 
the region of inhomogeneous phase(s), is different from 
that obtained from the MC simulations - compare Figs. 
[5] and Qjjb) - the A-line suggests a closed loop around 
the inhomogeneous phase(s). This seems to be inconsis- 
tent with the simulation results of Sec. IIII1 because the 
fluid of spherical clusters that we have identified is not 
periodically modulated. Thus while the RPA provides 
a reasonable approximation for the fluid structure and 
thermodynamics in regions of the phase diagram away 
from the A-line [32], it would appear to be unreliable in 
its direct vicinity. 



where —b(r) is the (unknown) bridge function. A simple 
approximation is to set b(r) = 0. This is the HNC closure 

c H Nc(r) = -Pv{r) + h(r) - ln(l + h(r)). (10) 

We solve Eqs. ® and (JTUJ) to obtain our approximation 
for g(r). Another approximation we consider is the PY 
closure |44j |: 

c PY {r) - (1 - exp[Pv{r)]){l + h(r)). (11) 

which is known to provide a good approximation for the 
structure and thermodynamics of a hard sphere fluid at 
low and moderate densities. 

Within the HNC approximation, the chemical poten- 
tial may be obtained from the following expression: 



B. Integral equation theories for the bulk fluid 
structure and thermodynamics 

The structure of a fluid may be characterized by the 
radial distribution function g(r) — l + h(r), where h(r) is 
the total correlation function, defined as the deviation of 
g(r) from the ideal vapor result [44j . The radial distribu- 
tion function is important for two main reasons: Firstly, 
it yields -via a Fourier transform-the static structure fac- 
tor: 



(3H = ln(pA 3 ) + p dr 



h(r) 



[fc(r)- C (r)]-c(rU,(12) 



S(k) = l + p / dr (g(r) - 1) exp(ik • r), 



(6) 



a quantity which, in principle, can be obtained in a scat- 
tering experiment. Secondly, g(r) can be used to calcu- 
late thermodynamic quantities such as the internal en- 
ergy and the pressure. For example, the pressure P can 
be obtained via the the virial equation [44j |: 



P = pk B T- £/dr ff (r)r^M 



(7) 



while, the Helmholtz free energy may be obtained by 
thermodynamic integration. In practice, however, since 
g(r) is only approximately known, the value obtained for 
the free energy is dependent upon the path of integration 

M 

A key equation used to calculate g(r) [or, equivalently 
h(r)] is the Ornstein-Zernike (OZ) equation [44j |: 

h(r) = c(r) + P J dr'cflr - r'\)h(r'), (8) 

where c(r) is the pair direct correlation function. To solve 
the OZ equation it must be supplemented by a closure 
relation. The exact closure reads [44| 

c(r) = -0v(r) + h{r) - ln(l + h{r)) - b(r), (9) 



where A is the (irrelevant) thermal de Broglie wavelength. 
Coexistence between two phases occurs if they have equal 
temperature, equal pressure and equal chemical poten- 
tial. We employed Eqs. and (fT2]) to calculate the 
fluid pressure and chemical potential to find the coex- 
isting densities within the HNC approximation. The re- 
sults are displayed in Fig. [T3J Unfortunately, only a small 
portion of the coexistence binodal (at low values of e _1 ) 
could be determined since for the HNC approximation 
there is a large region of the phase diagram in which 
there is no solution to the OZ equation. This region en- 
compases state points where we would expect (on the 
basis of the simulation) to find inhomogeneous phase(s) 
- see Fig. [121 We note that it is not unknown for the 
HNC approximation to exhibit regions of no solutions 
[59l IdOI. Kill lo2l|. Indeed, one can obtain the HNC closure 
(flTJ|) as the Euler-Lagrange equation from an associated 
free energy functional and it can be shown that the lack 
of HNC solutions is due to the loss of convexity in this 
free energy functional [U [12] ■ 

We have also applied a modified HNC (MHNC) closure 
in which we invoke a simple approximation for the bridge 
function in Eq. ©, chosen to be the PY expression for 
the hard-sphere fluid bridge function at the same density 
p. We found that this MHNC approximation failed to 
converge for the same state points as the HNC approx- 
imation. We did not attempt to implement any other 
more sophisticated MHNC approximations (58j . 

Applying the PY closure (TT|) to the OZ equation © 
gives some (perhaps) surprising results: As with the 
HNC, there are regions of the phase diagram where we 
were unable to obtain a solution. However, their location 
is different to that found for the HNC - compare Figs. [12] 
and [T31 The theory also seems to succeed in capturing 
some key aspects of the cluster transition. Specifically, 
if for some constant (T 1 < 0.48, one calculates the PY 
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FIG. 12: The phase diagram in the e _1 versus p plane ob- 
tained from the HNC integral equation theory for z\ = 2, 
Z2 = 1 and A = 0.55. The shaded region denotes state points 
for which no solution could be found. Outside this region, at 
low and moderate densities, the HNC theory provides a good 
approximation for g(r), when compared with MC data. We 
find that at low values of e' 1 the HNC predicts liquid-vapor 
coexistence. 



approximation for g(r) along a path of increasing den- 
sity, starting from very low densities, one finds at some p 
that a discontinuous jump occurs in the solution for g(r)— 
see for example the results in Figs. Q3] and [THl The lo- 
cus of points in the phase diagram at which g(r) 'jumps' 
in this way is displayed in Fig. [13] Note also, that if 
one traverses the reverse path, one can follow the high 
density branch of solutions to the PY displaying the clus- 
tering to very low densities. The low density solutions on 
this branch are particularly striking. The location in the 
phase diagram of the 'jump' line is close to where the 
transition to a cluster phase occurs, as determined from 
our MC simulations. 

In Fig. [TBlwe display \h(r)\ obtained from the PY the- 
ory for the state point e _1 = 0.41 and pa 3 — 0.220 (see 
also Fig. [T5|) . The logarithmic vertical axis allows inspec- 
tion of the long wavelength oscillatory decay of this cor- 
relation function, arising from inter-cluster correlations. 
A precursor to this oscillatory behavior is visible in our 
simulation results (Fig. [5]). We note that the amplitude 
of the oscillations decays rather slowly indicating very 
long ranged correlations in the cluster phase. 

We have attempted to calculate the binodal of the 
vapor-cluster transition from PY theory. Initially we 
sought to perform thermodynamic integration of the 
pressure (as obtained from the virial [7]) to acquire the 
Helmholtz free energy [58| and thence the coexisting 
state points. However, we were unsuccessful in this en- 
deavor due to insufficient self consistency within the the- 
ory. As an alternative, more approximate, approach we 
have done the following: i) We assumed that the 'jump' 
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FIG. 13: The phase diagram in the e _1 versus p plane as 
obtained by the PY integral equation theory. The shaded 
region denotes state points for which no solution could be 
found. The solid line indicates the line in the phase diagram 
at which there is a discontinuous change in g(r) - compare 
with Figs. [14] and 1151 This line is close the vapor-cluster 
phase transition in our MC simulations. The dot-dashed line 
indicates the locus of state point points at which the pressure 
[obtained via Eq. |J7J] equals that in the vapor phase. 



line obtained from the PY theory (see Fig. [T3"]) is roughly 
the coexistence curve for the vapor phase - our MC sim- 
ulation data supports this assumption, ii) To find the 
coexisting cluster phase density, we sought the cluster 
phase state point having the same value of e _1 and pres- 
sure [obtained via Eq. (J7JI] as the vapor phase at densities 
just below the jump line. This approach yields the dot- 
dashed line in Fig. 1131 Similar results can be obtained 
from an alternative approach, which is to plot the pres- 
sure obtained from the virial equation ([7]) as a function 
of density for fixed e _1 < 0.477; this produces a curve 
exhibiting a van der Waals loop (as well as a discontinu- 
ity at the jump line). Ignoring the discontinuity, one can 
perform the equal areas construction in order to obtain 
coexisting state points. 

In Fig. [T7] we display the static structure factor S(k) 
obtained from the PY theory for two state points either 
side of the vapor cluster-phase transition line. For these 
state points the PY result for g(r) agrees well with our 
MC simulation result - see Fig. [T5] S(k) displays a peak 
at small wave vector < k c <C 2n/a for states both 
sides of the transition line. This peak is known to be a 
signature of clustering. However, in the cluster phase the 
peak height is several orders of magnitude greater than 
in the vapor phase. 

Fig. [TH] shows the peak (maximal) value of S(k) ob- 
tained from the PY theory as a function of the fluid 
density p, for a selection of values of e _1 . At low p, 
S(k) ~ 1 Vfc. As the density is increased, the peak that 
grows the fastest, is that at k — k c , i.e. the curves in 
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FIG. 14: The radial distribution function, g(r), obtained from 
simulation (solid lines) and the PY theory (dashed lines) for 
(T x — 0.44 at the particular selection of densities indicated. 
Note that at pa = 0.058, the PY theory predicts a dramatic 
jump discontinuity in g(r). An increase in the magnitude of 
the first and second maxima also occurs in the MC simula- 
tion results at a slightly higher density, but it is much less 
pronounced than that predicted by the PY theory and does 
not appear to be discontinuous for this value of e _1 . 
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FIG. 15: The radial distribution function, g(r), obtained from 
simulation (solid lines) and the PY theory (dashed lines) for 
(T x — 0.41 at the particular selection of densities indicated. 
Note that in both the simulations and theory, a jump discon- 
tinuity occurs in g(r) between pa 3 — 0.024 and pa 3 — 0.109. 
Within the PY theory, this jump occurs at pa 3 = 0.04. 



Fig. [18] represent the dependence of S(k c ) on p in this 
regime. For values of e _1 < 0.48, a jump occurs in S(k c ) 
at the density of the vapor-cluster transition. Underly- 
ing this jump in the peak height is a wholesale discon- 
tinuous change in the entire form of S(k). As we have 
previously shown [3^], the value of the structure factor 
at k = is roughly proportional to its value at k = k c , 
and consequently, S(k = 0) also jumps at the transi- 
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FIG. 16: The modulus of the correlation function h(r) ob- 
tained from the PY theory for e" 1 = 0.41 and pa 3 = 0.220. 




FIG. 17: The static structure factor, S(k), obtained for the 
PY theory for the state points A = 0.55, e" 1 = 0.41 and for 
densities (a) pa 3 = 0.024, (b) pa 3 = 0.220. For both these 
state points the agreement between the PY results and the 
MC simulations for g(r) is good - see Fig. fTS] The inset to 
(b), shows S(k) on a logarithmic scale. 



tion point. Given, however, that 5(0) = pksT\T, where 
Xt is the isothermal compressibility, this in turn implies 
that a jump discontinuity occurs in xt at the transition 
density, in accordance with the simulation results. 

As the density is increased past the value at which S(k) 
jumps as a whole, the value of S(k c ) 7 ie. the peak height 
decreases. Eventually, at some higher density (the pre- 
cise value depending on the particular value of e _1 ), the 
maximal value of S(k) no longer occurs at k = k c =/= 0, 
but instead at k ~ 27r/cr, the next maximum in S(k). 
This next maximum originates from the hard-sphere cor- 
relations in the fluid, and the crossover in maximal value 
from one peak in S(k) to the other, accounts for the 
higher density discontinuity in the gradient of the curves 
displayed in Fig. 1181 For example, for the case when 



0.44, this occurs at pa 



0.6. 



In Fig. [TO] we display the direct pair correlation func- 
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FIG. 18: The behavior of the maximal value of the static 
structure factor S(k) as a function of p, obtained within the 
PY theory for a number of different values of e _1 , as described 
in the text. 
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FIG. 19: The direct pair correlation function c(r) obtained 
from the PY theory for the parameters A = 0.55 and e _1 = 
0.41 and density pa A = 0.220. The dashed line is a plot of 
—/3v(r). For the corresponding g(r), see Fig. 1151 and for the 
corresponding S(k), see Fig. 1171 

tion for the state point A — 0.55, e _1 = 0.41 and den- 
sity pa 3 — 0.220. For this state point the PY approx- 
imation for g(r) is in good agreement with the results 
from our MC simulations - see Fig. [151 Note that for 
1 < r/cr < 1.3, c(r) ^> —/3v(r). We believe it is for this 
reason that the MSA closure to the OZ equation: 

h(r) = — 1, r < er, 

c(r) = -0v(r), r > a, (13) 

completely fails to describe the fluid structure in the re- 
gion of the phase diagram near the vapor-cluster phase 
transition - the MSA approximation forces c(r) to be 
much smaller than it is in reality for 1 < r/a < 1.3. We 



believe that it is also for this reason that the RPA closure, 
Eq. (0, is unreliable in the vicinity of the vapor-cluster 
transition. 



V. DISCUSSION 

In summary, we have studied the phase behavior of 
a system of particles which interact via a mermaid po- 
tential, as a function of the strength of the short ranged 
attraction, for a fixed strength of the long ranged re- 
pulsion. The effect of the repulsion is to substitute the 
liquid- vapor critical point and a portion of the associated 
liquid-vapor transition line, with two first order phase 
transitions: one from the vapor to a fluid of spherical liq- 
uidlike clusters, and the other from the liquid to a fluid of 
spherical voids. Each of these phases is highly compress- 
ible compared to the respective homogeneous phase with 
which it coexists. At low temperature, the two transition 
lines intersect one another and the vapor-liquid transition 
line at a triple point, while at high temperatures, they 
appear to terminate at distinct critical points. Although 
SCOZA, HRT and most of the standard integral equa- 
tion theories are unable to describe the new transitions, 
somewhat surprisingly the Percus Yevick approximation 
does succeed in capturing the vapor-cluster phase transi- 
tion, as well as key aspects of the structure of the cluster 
fluid. 

Owing to the high computational cost of the current 
study, we have obtained the phase behaviour only for 
a single value of the strength of the repulsion namely 
A = 0.55. Consequently, it is difficult to comment au- 
thoritatively on the topology of the phase diagram in 
the full space of (j,, e and A. Nevertheless, it is wor- 
thy of note that the phase diagram of Fig. [5] is largely 
consistent with a cut at constant A through a phase di- 
agram originally obtained by Barbosa in a mean field 
study of an Ising model having isotropic competing in- 
teractions [63J. Fig. [2TjT a) recasts Barbosa's phase di- 
agram in terms appropriate for the present fluid model. 
One sees from the figure that as A is increased from zero, 
the liquid-vapor critical point shifts to lower e _1 in ac- 
cordance with SCOZA and HRT studies [3l[ of mermaid 
potentials. At some value of A = Al - corresponding 
to a Lifshitz point- inhomogeneous phases start to occur 
and two sheets (or wings) of first order phase transitions 
emerge from a triple line to delineate a region of inho- 
mogeneous states. This region is capped off at high e _1 
by a surface of continuous transitions. A cut through the 
phase diagram at constant A > Al has the form shown 
in Fig. I20f b). which is qualitatively very similar to our 
Fig. [SJ except that we have not, as yet, been able to 
identify a line of continuous transitions. 

A phase diagram of the form Fig. l2"0Ta) might shed 
light on the observation by Pini et al. of an anomalously 
large region of high compressibility around the critical 
point [3l|. If the strength of the repulsion is such that 
A < Al, then the system will exhibit a liquid- vapor crit- 
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FIG. 20: (a) A schematic representation of a possible form 
of the full phase diagram for the model of Eq. (TIJ, adapted 
from ref. [631 an d described in the text, (b) A cut at constant 
A> Al through the phase diagram of (a). 



ical point, but will be located (in the full phase diagram) 
close to the Lifshitz point. As we have seen in the present 
work, the inhomogeneous states that form for A > Al 
are highly compressible compared to the vapor or the liq- 
uid, and even though they aren't stable for A < Al, they 
would be expected to have an influence on the free energy 
landscape. This could account for the anomalously large 
region of high compressibility. 

As described in Sec. lIII Dl when traversing the inhomo- 
geneous region separating the low and high density tran- 
sitions, a variety of finite-size limited (spanning) struc- 
tures were observed depending on the value of e and 



the system size. Evidence was found for cylinders and 
slabs, as well as cylindrical voids. It is interesting to 
note that a similar sequence of inhomogeneous structures 
has also been reported in simulations of the sub-critical 
Lennard-Jones fluid using periodic boundary conditions 
[64j . There the observed structures occurred for values 
of the density intermediate between the coexisting stable 
vapor and liquid phases, and were metastable with re- 
spect to these phases. It therefore seems that the effect 
of adding a long ranged repulsion to a purely attractive 
fluid, is to stabilize these structures with respect to the 
homogeneous phases. We further note that while the 
clusters that occur at the low and high density transi- 
tions have a length scale smaller than our system size (at 
least for larger e _1 ), and are therefore expected to per- 
sist in the thermodynamic limit, our system sizes are too 
limited to make definitive statements regarding the sit- 
uation at moderate densities. Here it seems more likely 
that modulated structures occur, the wavelength of which 
exceeds our linear system sizes. 

With regard to the prospects for extensions to the 
present work, we were not able to satisfactorily ad- 
dress the question of the nature of the putative criti- 
cal points which terminate the vapor-cluster and liquid- 
cluster transitions. To do so will the require the identifi- 
cation of the appropriate order parameter for the transi- 
tions, and a reformulation of finite-size scaling method- 
ologies to take account of the large cluster length scale 
|65j |. Once this is achieved, the question as to whether 
a line of continuous transitions really does link the 
(tri)critical points (cf. Fig. I2tjfb)) could be tackled. As 
for future integral equation studies, another closure ap- 
proximation for the OZ equation, in common employ, is 
the Rogers- Young (RY) closure [58|, This interpo- 

lates between the PY and HNC theories in such a way 
as to enforce consistency between the virial and the com- 
pressibility routes to the thermodynamics. Since the PY 
theory is able to describe the vapor cluster-phase transi- 
tion, in contrast to the HNC, it would be interesting to 
see how the RY closure fares in this respect. 
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